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1 Basics 

Modern computer algebra systems use heuristics and algorithms for the fast com- 
putation with mathematical formulas. 

General purpose computer algebra systems like Axiom [O], Derive pD |, Mac- 



syma [lTj , Maple || , Mathematica [ 23j or Reduce ]TI) are for example great with 



integrals. Even a small system like Derive computes all explicitly given integrals 
of Bronshtein and Semedyayev's integral table |4j]. But how do such computations 
work? 

To begin with I would like to give examples of important mathematical con- 
cepts and methods that are available in general purpose computer algebra sys- 
tems. For demonstration purposes I use the Maple V.5 system. 

1.1 Linear Algebra 

One of the main topics of any computer algebra system is linear algebra. Linear 
algebra algorithms are used throughout Mathematics; we will see examples in 
connection with orthogonal polynomials and special functions later. 

With Maple, we can compute the solution of a linear system of equations: 

> solve ({x+2*a*y+3*z=4 , 5*x+6*y+7*z=8 , 9*x+10*y+l l*z=12} , 

> {x,y,z}); 

3 -1 

i z = 2> x = ~2> y = °) 

even if parameters are involved. For this purpose, Maple uses a Gauss type 
algorithm. 

Note, that the above system is linear only if a is considered constant. If we 
consider a as a variable, then a nonlinear system has to be solved: 
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> solve ({x+2*a*y+3*z=4 , 5*x+6*y+7*z=8 , 9*x+10*y+ll*z=12} , 

> {a,x,y,z}); 

3 —1 
{z = ^, x = —,y = 0, a = a},{ 

1 3 11 
z = -^V + ^, a = 1, x = --y - -, y = y} . 

In a forthcoming section, we give more details on nonlinear systems of equations. 
Maple has a large linear algebra library: 

> with(linalg) ; 

Warning, new definition for norm 
Warning, new definition for trace 

[BlockDiagonal, GramSchmidt, JordanBlock, LUdecomp, QRdecomp, 
Wronskian, addcol, addrow, adj, adjoint, angle, augment, backsub, 
band, basis, bezout, blockmatrix , charmat, charpoly, cholesky, col, 
coldim, colspace, colspan, companion, concat, cond, copyinto, 
crossprod, curl, definite, delcols, delrows, det, diag, diverge, 
dotprod, eigenvals, eigenvalues, eigenvectors, eigenvects, 
entermatrix, equal, exponential, extend, ffgausselim, fibonacci, 
forwardsub, frobenius, gausselim, gaussjord, geneqns, genmatrix, 
grad, hadamard, hermite, hessian, hilbert, htranspose, ihermite, 
indexfunc, innerprod, intbasis, inverse, ismith, issimilar, iszero, 
jacobian, jordan, kernel, laplacian, leastsqrs, linsolve, matadd, 
matrix, minor, minpoly, mulcol, mulrow, multiply, norm, normalize, 
nullspace, orthog, permanent, pivot, potential, randmatrix, 
randvector, rank, ratform, row, rowdim, rowspace, rowspan, rref, 
scalarmul, singularvals, smith, stackmatrix, submatrix, subvector, 
sumbasis, swapcol, swaprow, Sylvester, toeplitz, trace, transpose, 
vandermonde, vecpotent, vectdim, vector, wronskian] 

You can see which procedures are available now. As an example, we compute the 

determinant of the matrix 

/ 1 2a 
4 5 
\ 7 8 

by 

> det([[l,2*a,3] , [5,6,7] , [9,10,11]]); 

-16 + 16a 

and the eigenvalues and eigenvectors for a = 1: 

> eigenvalues ( [[1,2,3] , [5,6,7] , [9,10,11]]) ; 

o, 9 + vTos, 9 - VToE 
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> eigenvectors ( [[1,2,3] , [5,6,7] , [9, 10, 11]] ) ; 

[0, 1, {[1, -2, 1]}],[9 + Vl05, 1, 

{ 



— -- >/l05, 1, -- + - VT05 
2 2 ''22 



9 - vTos, 1, 



{ 



— + - Vl05, 1, -- - - VT05 
2 2 ' ' 2 2 



}]• 



}] 
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Later we will show how important an efficient implementation of linear algebra 
can be. 

1.2 Polynomial Artithmetic 

A second major topic of computer algebra is polynomial arithmetic. 
P defines a polynomial 

> P:=(l-x)*sum(x~k,k=0. .9) ; 

P ■= (1 - x) (1 + x + x 2 + x 3 + x 4 + x 5 + x 6 + x 7 + x 8 + x 9 ) 

which is expanded by 

> expand (P) ; 

Q is a second polynomial 

> Q:=l-x~5; 

Q := 1 - x 5 

and normal cancels common factors of the ratio P/Q: 

> normal (P/Q) ; 

x 5 + l 

by an application of the Euclidean algorithm. 

A highlight of computer algebra is rational factorization since the underly- 
ing algorithms (factorization modulo a prime and Hensel lifting, or the triple L 
algorithm, see e.g. [§]) are not suitable for hand computations. 

For example, the polynomial P can be factored over Q by the command 

> factor(P) ; 

-(-1 + x) (x + 1) (x 4 + x 3 + x 2 + x + 1) (x 4 - x 3 + x 2 - X + 1) 
and the following is a rational factorization of 1 — x 105 : 

> factor(l-x"105) ; 



-(-1 + x) (x 6 + x 5 + x 4 + x 3 + x 2 + x + 1) (x 4 + x 3 + x 2 + X + 1)(1 - X + x 

_ ™6 i 7 _ 8 , 10 _ 11 , _12 _ „13 , 14 _ 16 , 17 _ 18 , 19 

- x 23 + x 24 )(x 2 + x + 1) (1 - x + x 3 - x 4 + x 6 - x 8 + x 9 - x 11 + x 12 ) 
(1 - x + x 3 - x 4 + x 5 - x 7 + x 8 ) (1 + x + x 2 - x 5 - x 6 - 2 x 7 - x 8 - x 9 

_ 24 i 12 , 13 , 14 , 16 , 17 , ~15 , 48 _ 20 _ 22 _ 26 _ 28 

iAj \ iAj \ JU \ *Xj \ JU \ iAj \ *Xj \ JU iAj iAj iAj JU 
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+ x 31 + ^32 + ^33 + ^34 + ^35 + ^36 _ ^39 _ ^40 _ 2 % Al _ ^42 _ ^43 + ^46 

+ x 47 ^ 
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Note that 105 is the smallest exponent such that the rational factorization of 
1 — x n contains coefficients different from or ±1. 
Next, we define a multivariate polynomial 

> Product (x~(2*k-l)-y~k/k~2,k=l. .7) ; 

k=l 

whose expanded form is a huge expression: 

> term :=expand (product (x~ (2*k-l)-y~k/k~2 ,k=l . . 7) ) ; 

term := _1 x ™ 7 _ 1 40 5 + J_ ^ 12 + _i 2 l_ 3.29 11 
49 y 25 y 1225 y 176400 y 

' ,16 18 1 ™42„.4 , 31 10 1 18 17 



x 2/ ~ 77 x 2/ + 7^777 x 2/ ~ 777777 x V 



44100 16 28224 a 28224 

16969 33 9 49 _ 181 20 16 _ 71 22 15 

1587600 y 1587600 y 235200 y 

+ X 9 y 22 - - X 44 y 3 + x 35 8 _ ^24 14 

705600 y 9 y 176400 y 45360 y 

396900 y 3600 y 3175200 y 6350400 y 
2069 y i2 + _J_ x w y i9 + _1 3 ^ ^ 18 1 ^ 4 y25 



564480 26880 181440 6350400 

I x 46 t/ 2 + 89 x 39 t/ 6 - 17147 x 30 v 11 + 559 X 19 7/ 17 
4 y 1600 y 3175200 y 3628800 y 

1 x 6 24 , i3_ 41 5 _ 18733 32 10 , 167 21 16 

2822400 y 144 y 1587600 y 453600 y 

13 3.8 23 _ 21 ^34 9 , 26 29 x 23 15 

6350400 y 1600 y 5080320 y 

89 3.10 22 , 1091 25 14 _ 1 12 21 

25401600 y 1270080 y 90720 y 

J4„.20 i 1 „„.27 „.„.48 i 1 ,4 43 



x ± - ± y™ + xy" -yx w + -y x 



12700800 * 25401600 * 9 

_ w 8 ^,36 , 401 13 ^27 2 j^_ 19 16 , 1 26 3 

3600 y 313600 y 635040 y 6350400 y 

+ I u 3 45 _ J_ 7 38 , 181 12 29 113 ,.18 18 

+ 4^^ 64 y + 129600 V 2822400 y 

+ I y 2 ^ _ J_ 6 40 , J_ 11 31 421 y 17 20 

2822400 y 36 y 900 y 6350400 y 

+ \ y24 3.7 , J_ 10 3,33 L_ yl6 3,22 , 1 23 x 9 

1587600 576 20736 1016064 y 

_ 1 x ™ y 6 + J_ x 2 5 13 _ _J_ 15 ^24 + _L_ 22 11 

36 1764 y 14400 705600 

+ _^ 2 / 2 V 3 y 28 

518400 y 25401600 y 
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which is a polynomial of degree 49 w.r.t. x and of degree 28 w.r.t. y. 

It is beautiful (and will turn out to be essential in the sequel) that computer 
algebra systems have no problems to factorize such expressions over the rationals 
in reasonable time: 

> factor (term) ; 



1.3 Polynomial Systems 

We come back to the problem of nonlinear systems of equations. Whereas in the 
linear case, Gauss elimination works, Buchberger's algorithm is an extension to 
the multivariate case. It constitutes — given a certain term order — an elimination 
scheme to find a normal form for a given polynomial system, which can be used 
to find the general solution of a nonlinear system. 
We consider the following system of equations: 

> LIST:={9*B*A+4*d-6*c*d=0, 

> -9*a*b+9*B*A=0, 

> -18*B*A+12*d-12*c*d+4*d~2=0, 

> 6*b*d-36*a*b+2*d+6*a*d=0, 

> -4*d~2+12*b*d-36*a*b+12*a*d=0, 

> -8*C-9+9*B+9*A-4*d+12*c=0, 

> -8*d-7+12*c+27*A+27*B-24*C=0, 

> 8-3*a-3*b-32*C+27*A+27*B=0, 

> 6-16*C+18*A+18*B-12*a+4*d-12*b=0 , 

> 4-12*a-12*b+8*d=0, 

> -C-2+3*c=0}; 

LIST := {9BA + 4d-6cd, -9ab + 9BA, 

-18B A + 12 d- 12 cd + 4d 2 , 6bd- 36 ab + 2 d + 6 ad, 
-4 d 2 + 12 b d - 36 a b + 12 a d, 4 - 12 a - 12 b + 8 d, 
-C-2 + 3c, -8C-9 + 9B + 9A-4d+12c, 
-8d-7 + 12c + 27 A + 27 B -24C, 
8-3a-36-32C + 27 A + 27 B, 
6 - 16 C + 18 A + 18 B - 12 a + 4 d - 12 b} . 
The solve command gives the general solution: 

> solve(LIST,{A,B,C,a,b,c,d}) ; 




— (x - y) {-y 2 + 4 x 3 ) (-y 3 + 9 x 5 ) (-y 5 + 25 x 9 ) {-y G + 36 a; 11 ) 
(49 x 13 -y 7 ) (-y 4 + lQx 7 ) 




C = ^ + -d,d = d}, 

1 , n 1 1 j 
- a, C = - + - a, 



5 




C = - + -d, d = d\,\c=- + -d, b=-d+-,B 
2 2' I,x 6 6' 3 3' 



1 , 1 A 1 j 

— cl ~\ — , A = - d, 
3 3' 3 ' 



a = - d, C = — I — d, d — d\ . 



In an application, we will meet this example later again. 
1.4 Differentiation and Integration 

Differentiation is done using the differentiation rules. This is an easy task. For 
our example function 

> input : =exp(x-x~2) *sin(x~6-l) ; 



obviously the product rule is used: 

> derivative :=diff (input ,x) ; 

derivative := (1 — 2 x) e^ x ~ x ^ sin(a; 6 — 1) + 6 e^~ x ^ cos(a; 6 — 1) x 5 . 
Integration is much more difficult, and the different systems have different ap- 
proaches: Whereas Derive uses a good collection of heuristics which enable the 
system to compute all explicitly given integrals of Bronshtein and Semedyayev's 
integral table [|J], as already mentioned, Maple uses an algorithmic approach. 

In the sixties Risch developed an algorithm to compute an elementary an- 
tiderivative whenever one exists. If no such antiderivative exists, his algorithm 
returns this information. Here elementary means that both integrand and an- 
tiderivative are rationally composed of exponentials and logarithms (see e.g. ||). 
Adjoining the complex unit i (denoted in Maple by /), trigonometric functions 
can be treated as well. 

We integrate the derivative above. This takes a little longer: 

> integral : =int (derivative , x) ; 



Since i is adjoined, the resulting function looks not very familiar although it is 
algebraically equal to our input function. In this particular case, we can con- 
vert both functions to the same normal form by first converting exponentials to 
trigonometries and applying then rational factorization: 
> factor(convert(integral,trig)) ; 



input := e 



sm(x — 1) 



integral : 




sin((x - 1) (x + 1) (x 2 + x + 1) (x 2 - x + 1)) 
(— cosh(x (x — 1)) + sinh(x (x — 1))) 
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> f actor (convert (input, trig) ) ; 

-sin((x - 1) (x + 1) (x 2 + x + 1) (x 2 - x + 1)) 
(— cosh(x (x — 1)) + sinh(x (x — 1))) 
Note, however, that one can prove that for general transcendental expressions a 
normal form does not exist. 



1.5 Differential Equations 

In engineering and in natural sciences the symbolic and numeric solution of differ- 
ential equations is rather important. We enter an ordinary differential equation: 

> DE:=diff (y(x) ,x)=l+y(x) "2; 

DE :=^-y(x) = l+y(x) 2 . 

After loading the DEtools package, we can use the procedure df ieldplot to plot 
a direction field of the differential equation: 

> with(DEtools) : 

> df ieldplot (DE,y(x) ,x=-5. .5,y=-5. .5) ; 




Given an initial value, the command DEplot plots a numeric solution by a Runge- 
Kutta type approach: 

> DEplot ({DE},{y(x)},x=-l. .1, [ [y(0)=0] ] ) ; 
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Using a combination of heuristic and algorithmic techniques, Maple can solve 
many ordinary differential equations explicitly. Our initial value problem has the 
solution: 

> dsolve({DE,y(0)=0},y(x)); 

y(x) = tan(x) . 

As another example, we consider a linear differential equation of second order. 

> DE:=diff (y(x) ,x$2)-y(x)=sin(x)*x; 

DE : = (— y(x)) - y(x) = sin{x) x 

with explicit solution 

> dsolve(DE,y(x)) ; 

III l III 

y( x ) = y~2 Z+ 2^ X C ° S ^ + 4 Sin ^' ) X ^ + 2 *~ 2 X " 2 ^ ei ~ X) C ° S( ^ 
I III I 

— - x e^~ x ^ sin(:r))sinh(:r) + (— - (— - x + -) e x cos(x) — - sin(x) x e x 

+ - (— - x — -) e*"^ cos(x) — - xe^~ x ^ sin(a;))cosh(a;) + _C1 sinh(:r) 

+ _C2 cosh(x) . 
A plot based on a numerical computation is given by 

> DEplot({DE},{y(x)},x=-5. . 5 , [ [y (0) =0 ,D(y) (0)=1] ] ) ; 
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The corresponding initial value problem has the explicit solution 
> solution:=dsolve({DE,y(0)=0,D(y) (0)=l},y(x)) ; 



solution := y(x) 



'- ( — x H — ) e x cos(x) H — sin(rr) x e x 



2 v 2 



1 

4 



+ - (— - x — — ) e 1 - ^ cos(a;) --ie' ^ sin(x))sinh(a;) + ( 

(-^ z + g ) e * cos ( x ) ~ \ sin ( x ) x e * + 2 ^ X ~ ^ e(_:E) COs( ^ 

— -ie' _I ' sin(x))cosh(x) + sinh(x) + - cosh(x) 

which can be simplified to 

> simplify (convert (rhs (solution) ,trig)) ; 

sinh(x) H — cosh(x) cos(x) sin (a;) x . 



1.6 Formal Power Series and Differential Equations 

Next, we consider the opposite problem to generate differential equations from 
expressions. This will lead us also to the generation of power series of hypergeo- 
metric type. 



After loading the FPS package [10| 
> with(share) : with(FPS) : 
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See ?share and ?share , contents for information about the share library 
Share Library: FPS 
Author: Gruntz, Dominik. 

Description: FPS function attempts to find a formal power 
series expansion for a function in terms of a formula for the 
coefficients 

we can, e.g., compute the formal power series of the square of the inverse tangent 
function: 

> FPS(arcsin(x)~2,x) ; 

~ ( k[) 2 A k x (2k+2) 

^(* + l)(l + 2*) r U 



The algorithm behind this procedure is the following ([14|, ||10|| ): 

In the first step, by linear algebra techniques, a homogeneous linear differential 
equation with polynomial coefficients is sought for the given expression 

> DE:=SimpleDE(arcsin(x) ~2,x,F) ; 

DE := (x -l)(x+ 1) F(x)) + F(x)) + 3 x F(x)) = . 

We call such a differential equation as well as the corresponding function holo- 
nomic. Next, substituting the series 

oo 
k=0 

in this differential equation and equating coefficients yields the holonomic recur- 
rence equation for a^: 

> RE:=SimpleRE(arcsin(x)~2,x,a) ; 

RE :— —(k + l)(a(Jfe + 3) k 2 - k 2 a(fc + 1) + 5 a(Jfe + 3) k - 2 k a(fc + 1) 
-a(fc + 1) + 6a(A; + 3)) =0 
which can be put in factored form 

> map(f actor, collect (lhs (RE) ,a))=0; 

-(k + l)(k + 2) (k + 3) a(k + 3) + (k + l) 3 &{k + 1) = . (2) 
Notice that the resulting recurrence equation gives a k +2 as a rational multiple 
of a*;. If Ak + i is a rational multiple of A k then it is called a hyper geometric 
term. From (0), can be easily computed using two initial values. This finally 
generates the explicit series representation ([]]). Note, however, that for an explicit 
representation the above factorization is necessary; see (0). 

By solving the differential equation for F(x) = arcsin 2 (a;) with two initial 
values, we would like to reconstruct the input: 

> solution:=dsolve({DE,F(0)=0,D(F) (0)=0, (D@@2) (F) (0)=2},F(x)) ; 
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solution := 

F(x) = - 7T 2 + / 7T ln(x + y/(x ~ 1) (x + 1)) - hl(x + ^(x - 1) (x + l)) 2 

> convert (arcsin(x) "2 , In) ; 

-ln(Vl -x 2 + Ixf . 
As before, we see that transcendental functions come in quite different disguises. 

It turns out that sum and product of two holonomic functions are again 
holonomic, and the corresponding holonomic (differential or recurrence) equations 
can be constructed from the given holonomic equations by linear algebra ( |lj , , 

As an example, we consider both the sum and the product of the functions 
f(x) = arcsinx and g(x) = e x . Here are their holonomic equations: 

> DEI : =SimpleDE(arcsin(x) ,x,F) ; 

DEI := {x - 1) (x + 1) F(x)) + F(ar)) x = 

ox ox 

> DE2:=SimpleDE(exp(x) ,x,F) ; 

.9 

DE2 := (—F(x))-F(x) = . 
ox 

From these, we can compute the holonomic equations that are valid for f(x)+g(x) 



and f(x) ■ g(x). For this purpose, we load the gfun package |2T 

> with (gfun) ; 

[Laplace, algebraicsubs, algeqtodiffeq, algeqto series, algfuntoalgeq, hovel, 
cauchyproduct, diffeq * diffeq, diffeq + diffeq, diffeqtorec, guesseqn, 
guessgf, hadamardproduct, holexprtodiffeq, invborel, listtoalgeq, 
listtodiffeq, listtohypergeom, listtolist, listtoratpoly , listtorec, 
listtoseries, listtoseries/ Laplace, listto series j egj , listtoseries / Igdegf , 
listtoseries / Igdogf , listtoseries / ogf , listtoseries / revegf , 
listtoseries / rev ogf , maxdegcoeff , maxdegeqn, maxordereqn, 
mindegcoeff , mindegeqn, minordereqn, optionsgf, poltodiffeq, 
poltorec, ratpolytocoeff , rec * rec, rec + rec, rectodiffeq, rectoproc, 
seriestoalgeq, seriestodiffeq, seriestohypergeom, seriestolist, 
seriestoratpoly , seriestorec, seriestoseries] 
The procedures 'dif f eq+dif f eq' and 'dif f eq*dif f eq c compute the differen- 
tial equations of sum and product, respectively: 

> 'diffeq+diffeq' (DEI ,DE2 , F (x) ) ; 

{(-x 3 -2x 2 +x- l)D(F)(x) + (-x 4 + 4x 2 ) (D (2) )(F)(x) 

+ (1 - 2x 2 + x 4 - x + x 3 ) (D( 3 ))(F)(x), (D( 2 ))(F)(0) = _C } 

> 'diffeq*diffeq f (DEI ,DE2 , F (x) ) ; 

(-1 + x 2 - x) F(x) + (x + 2 - 2x 2 ) D(F)(x) + (-1 + x 2 ) (D (2) )(F)(x) 
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which we could also have obtained using SimpleDE directly:^ 
> SimpleDE(arcsin(x)+exp(x) ,x,F) ; 

d 3 d 
(x - 1) (x + 1) (x - 1 + x 2 ) (— F(x)) + ix - x 3 - 2 x 2 - 1) (— F(x)) 

cte d ox 



d 2 

-x 2 (x-2) (x + 2)(— F(x)) = 

> SimpleDE(arcsin(x)*exp(x) ,x,F) ; 

(x - 1) (x + 1) F(x)) + (x + 2 - 2 x 2 ) F(x)) 
ox z ox 

+ (-1 +x 2 -x)F(x) = . 

SimpleDE can also generate differential equations for some special functions, e.g., 

for the Bessel functions J n {x): 

> DE:=SimpleDE(BesselJ(n,x) ,x,F) ; 

d 2 d 
DE := (— F(x)) x 2 - (n - x) (n + x) Fix) + (— Fix)) x = . 
ox 2 ox 
Maple can solve this differential equation easily: 

> dsolve(DE,F(x)) ; 

F(x) = -CI BesselY(n, x) + _C^BesselJ(n, x) . 
Even the more complicated differential equation of the product 

> DE:=SimpleDE(BesselJ(n,x)*exp(x) ,x,F) ; 

d d 2 

DE := (2x 2 -n 2 - x)F(x) - (-1 + 2 x) x (— F(x)) + {—F{x))x 2 = 

ox ox 2 

can be treated by Maple 

> dsolve(DE,F(x)) ; 

F(x) = _Ci BesselJ(n, x) e x + _C2 BesselY(n, x) e x , 
but for the differential equation 

> DE:=SimpleDE(BesselJ(n,x)+exp(x) ,x,F) ; 

DE := (2 x 4 - 3 x 2 n 2 + x 3 - 3 n 2 x - x 2 + n 4 - n 2 ) F(x) 

d 

+ i-n 4 + 3 x 2 + n 2 + x 3 - 2 x 4 + 3 x 2 n 2 ) i— Fix)) 

ox 

d 2 

- (-2 x 3 + n 2 x + x 2 - 3 n 2 + 2 x) x i— - Fix)) 

ox 2 

d 3 

+ x 2 (-2x 2 + n 2 -x)(— F(x)) = 0, 

ox 6 



1 Note that SimpleDE uses a slightly different approach (also based on linear algebra) 
that sometimes can find differential equations of lower order than 'dif f eq+dif f eq' and 
'diffeq*diffeq'. 



12 



Maple fails: 

> dsolve(DE,F(x)) ; 



F(x) = _C1 e x + e^DESol 



x {2 x + !){§- _Y(x)) (-2 a; 4 + a; 2 n 2 - x 3 ) _Y(x))' 



— 2 a; 2 + n 2 — x (—2 x 2 + n 2 — re) 2 

{_Y(x)}J 

although Maple was able to find the exponential summand (and hence reduced 
the order by one). 

Nevertheless, it is not astonishing that Maple cannot find all such solutions 
since for this type of nonelementary solutions no algorithms exist. 



2 Special Functions and Computer Algebra 

Power series of hypergeometric type — the example function arcsin 2 x as well as the 
Bessel functions are of this type, e.g. — are the most important special functions. 
The generalized hypergeometric series is given by 



p 1 q 



di a 2 ■ ■ ■ % 
bi b 2 ■■■ b q 



~h(bi) k -(h) k ---(b q ) k k\ X [} 



k 

where (a)k '■= Yl ( a +J — 1) = r(a + k)/T(a) denotes the Pochhammer-Symbol or 
j i 

shifted factorial. 

Ak is a hypergeometric term and fulfils the recurrence equation (k G N) 

(k + Oi) ■ ■ ■ (fc + Op) 
(fc + 6 1 )...(ib + 6,)(A ; + l) fc 

with the initial value 

A := 1 . 

In Maple the hypergeometric series is given as hypergeom(plist,qlist,x), 
where 

plist = [ai, a 2 , . . . , a p ] , 
qlist = [6i,6 2 , • • • ,b q ] . 

Here are some more hypergeometric examples: 

> F:=sqrt(x)*arcsin(sqrt(x))+sqrt(l-x) ; 



F := \[x arcsin {yfx) + yl 



x 
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> SUM:=FPS(F,x) ; 

SUM . y 4(-*)(2ifc)!^ 
S(*0 8 (2*-l) a 

> convert (SUM, hypergeom) ; 

hypergeom([— , — ], [-], x) 

> F : =- (sqrt (Pi) /2*sqrt (x) *erf (sqrt (x) ) * ( 1+1/2/x) +exp (-x) /2) ; 

F := -\ V5F Vxerf(Vi) (1 + ^) - ^ e<-> 

> SUM:=FPS(F,x) ; 

^k\ (2* + l)(-l + 2*) 

> convert (SUM, hypergeom) ; 

— KummerM(— , -, — x) 

With convert, one can convert series into hypergeometric notation; KummerM is 
another name for the confluent hypergeometric function iF\. 

2.1 Summation 

Whereas the FPS command converts expressions into series representations, the 
opposite question is to find explicit representations for sums. Note that the 



examples of the remaining paper are from ||15|| . 

The main interest lies in sums of hypergeometric terms. As an example, we 
ask: Why does Maple evaluate the sum 

> sum((-l)~k*binomial(n,k) ,k=a. .b) ; 

(6 + 1) binomial(n, 6 + 1) a (-l) a binomial(n, a) 

n n 
for arbitrary bounds a and b in simple form, but fails with 

> sum (binomial (n,k) ,k=a. .b) ; 

binomial(ra, a) hypergeom([l, — n + a], [1 + a], —1) 

— binomial (n, 6 + 1) hypergeom([l, — n + 6+1], [2 + 6], —1) ? 
On the other hand, for the special bounds a = and 6 = n, Maple is successful, 
again: 

> sum (binomial (n, k) ,k=0. .n) ; 

2 n . 

The reason for this behavior is that the first summand, (— l) k (?J, has a hyperge- 
ometric term antidifference (w.r.t. the variable k), and the second one, (?J, has 
not. The last sum is a definite sum with natural bounds, i.e., the sum can be 
considered as infinite sum (k = — 00..00), and the result, again, is a hypergeo- 
metric term (w.r.t. the variable n). We will see how we can find these types of 
results algorithmically. 
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Sk is called an antidifference of a*,, if 



Sk+l — Sk — a k , 

If such an antidifference is known, then summation is trivial since by telescoping 



b 



^ a k = - «&) + (Sb - «6-l) H 1" (Sa+1 - S a) = _ S a ■ 

k=a 

This is very similar to the integration case. 

The antidifference of the first summand is given by 

> sum((-l)~k*binomial(n,k) ,k) ; 

k (— l) k binomial(n, k) 

n 

We can increase the level of user information by the command 

> inf olevel [sum] : =3 : 

Let's try yo prove the statement 

(4jfc + l)(2Jfe)! 



^ k\A k {2k- 1) (ife + 1)! 



that was posed in SIAM Review 36, 1994, Problem 94-2 | 19| . We compute an 
antidifference 

> summand:=(-l)-(k+l)*(4*k+l)*(2*k) ! /(k ! *4~k* (2*k-l) *(k+l) !) : 

> sum ( summand, k) ; 

sum/indefnew: indefinite summation 

sum/extgosper : applying Gosper algorithm to a( k ) := 
(-l)~(k+l)*(4*k+l)*(2*k) !/k!/(4"k)/(2*k-l)/(k+l) ! 
sum/gospernew: a( k )/a( k -1):= 
-1/2* (4*k+l) / (4*k-3) / (k+1) * (2*k-3) 
sum/gospernew: Gosper's algorithm applicable 

4*k+l 
-2*k+3 
2*k+2 

sum/gospernew: degreebound:= 
sum/gospernew: solving equations to find f 
sum/gospernew: Gosper's algorithm successful 
sum/gospernew: f:= -1 

sum/indefnew: indefinite summation finished 

(k + 1) (-l)( fc +!) (2k)\ 

~ 2 k\A k (2k-l) (fc + 1)! ' 
with success. Taking the limit as n — > oo, one gets therefore 

> sum ( summand, k=l . . infinity) ; 

sum/inf inite : infinite summation 

1 . 



sum/gospernew: p 
sum/gospernew: q 
sum/gospernew: r 
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Moreover, from the user information we see that Gosper's algorithm is applied. 
If a k is a hypergeometric term, i.e., iff] 

then Gosper's algorithm decides whether or not the antidifference s k is a hyper- 
geometric term, and computes it in the affirmative case. 
In detail: Given[| 



a k+l _ bk 

a representation 



bk _ Pk+1 Qk+l 

is computed for which 



c k p k r k+ i ' 



b k) c k E Q[k] 



Pk,Qk,r k e 



gcd (q k , r k+j ) = 1 for all j G N . 

This can be done by a resultant computation || or by rational factorization ([0. 
[0)- 

The essential fact is then: f k , defined by 

r k 

Sk — — Jk-l 0, k 

Pk 

is rational, and the above gcd-condition yields even f k e Q[k}. f k satisfies the 
inhomogeneous recurrence equation 

p k = q k+1 f k - r k . 

After calculating the degree of f k , it is pure linear algebra to compute f k . The 
output of the procedure is either 

r k , 

Sk — — Jk-l a k 
Pk 

or the statement "There is no elementary (= hypergeometric term) antidiffer- 
ence". 

In the book fl5(| , many algorithms that are connected with Gosper's, are 
treated in detail and Maple implementations are given. 
After loading 1 hsum . mpl ' ,Q 
> read ( 'hsum. mpl ' ) ; 

Copyright 1998 Wolfram Koepf, Konrad — Zuse — Zentrum Berlin 



2 Q(k): rational functions over Q. 
3 <Q>[fc]: polynomials over Q. 
The packages 'hsum. mpl' and 'qsum.mpl' can be obtained from the URL 
www . imn . htwk-leipzig . de/ "koepf /research . html. 
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we can repeat the above calculation by the command 

> gosper((-l)~(k+l)*(4*k+l)*(2*k) !/(k! *4~k* (2*k-l) * (k+1) D ,k); 

_ (fc+l)(-l)( fc+1 ) (2k)\ 
k\4 k (2 k — 1) (k + T)\ ' 

The computation 

> gosper(l/k,k) ; 

Error, (in gosper) no hypergeometric term antidif f erence exists 

is not worthless at all: It proves that the harmonic numbers 



n 1 



k 

k=l 

do not constitute a hypergeometric term, corresponding to the fact that the 
logarithmic function 

r 1 

In x = / - dt 

cannot be written in terms of exponentials. 

Zeilberger's algorithms is an extension of Gosper's for definite sums. It gen- 
erates, e.g., the right-hand sides of the identities 

£(;)-*■ 

fc=0 v 7 

v ( n Y _ ( 2n ) ! 

by the commands 

> closedf orm (binomial (n,k) ,k,n) ; 

T 

> closedform(binomial(n,k)~2,k,n) ; 

(2n)\ 

(n\f ' 

Here are the details: If F(n, k) is a hypergeometric term w.r.t. n and k, i.e. 

F(n + l,k) , F(n,A;+l) (rfc . 
L , \ and 1 J €Qn,fc , 

F(n, k) F(n, k) 

then Zeilberger's algorithm generates a holonomic recurrence equation for 

s n := J^F(n, k) . 

kez 
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This is performed by starting with J = 1 and iterating if necessary: Set 

j 

a k := F(n, k) + ^ <7j(n) F(n + j, k) 
3=1 

with as yet undetermined variables a r Apply Gosper's algorithm to a k . In the 
last step, solve the linear system at the same time for the coefficients of f k and 
the variables Oj (j — 1, . . . , J). In the affirmative case, this yields 

G(n, k + 1) — G(n, k) = a k ■ 

Output: By summation one gets: 

j 

s n + ^2<?j(n) s n+j = . 

3=1 

We would like to point out that the most time consuming part of Zeilberger's 
algorithm is its last step which is to solve a linear system. This linear system, 
however, often has many variables, and its coefficients are polynomials or ratio- 
nal functions. Here, an efficient implementation of linear algebra is important. 
Furthermore, the resulting recurrence equation usually needs factored coefficients 
because otherwise the results look unnecessarily complicated. We will see such a 
situation soon. 

We give some examples: Each of the following series represents the Legendre 
polynomials: 
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Again, you see, that functions come in quite different disguises. How can we show 
that these systems define the same family of functions? Zeilberger's paradigm is 
to show that they satisfy the same recurrence equation, then it is sufficient to 
check a finite number of initial values. 

Here are the recurrence equations for the different sums: 

> P:='P': 

> sumrecurs ion (binomial (n,k) *binomial (-n-1 ,k) * ( (l-x)/2) ~k,k,P(n) ) ; 

(n + 2) P(n + 2) - (2 n + 3) x P(n + 1) + (n + 1) P(n) = 

> sumrecursion(l/2~n*binomial(n,k) ~2*(x-l) " (n-k)*(x+l) ~k,k,P(n)) ; 

(n + 2) P(n + 2) - (2 n + 3) x P(n + 1) + (n + 1) P(n) = 

> sumrecursion(l/2~n*(-l) ~k*binomial (n,k) * 

> binomial (2*n-2*k,n) *x~ (n-2*k) ,k,P(n)) ; 

(n + 2) P(n + 2) - (2 n + 3) x P(n + 1) + (n + 1) P(n) = 

> sumrecursion(x~n*hyperterm( [-n/2,-n/2+l/2] , [1] , l-l/x"2,k) ,k,P(n)) ; 

(n + 2) P(n + 2) - (2 n + 3) x P(n + 1) + (n + 1) P(n) = 
We omit the computation of the initial values. 

The Sumtohyper procedure of the hsum package is slightly more efficient than 
' convert /hypergeom' by converting a series into hypergeometric notation. 

> Sumtohyper (binomial (n,k) *binomial (-n-1 ,k) * ( (1-x) /2) ~k,k) ; 

11 

Hypergeom([n+ 1, -n], [1], - - -x) 

> Sumtohyper (l/2"n*binomial(n,k) ~2*(x-l) " (n-k)*(x+l) ~k,k) ; 

11 x + 1 

(-x- -) n Hypergeom([-n, -n], [1], ^— y) 

> Sumtohyper (l/2"n*(-l) "k*binomial(n,k)* 

> binomial (2*n-2*k,n) *x~ (n-2*k) ,k) ; 

2 ( ~ n * ) binomial(2 n, n) x n Hypergeom([ — n-\ — , — n], [— n H — ], — ) . 

The above computations show that all the given representations of the Legendre 
polynomials agree. 

To give a more advanced example of an application of Sumtohyper, we com- 
pute the hypergeometric representation of the difference P n+ i(x) — P n (x) of suc- 
cessive Legendre polynomials: 

> legendreterm: =binomial (n,k) *binomial (-n-1 ,k) * ( (l-x)/2) ~k; 

legendreterm := binomial(n, k) binomial(— n — 1, k) (- — - x) k 

> Sumtohyper (subs (n=n+l, legendreterm) -legendreterm, k) ; 

(x + xn — 1 — n) Hypergeom( [— n, n + 2], [2] , - — - x) . 
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We give more examples of how Zeilberger's algorithm can be applied in rather 
different situations. 

The following recurrence equation of the Apery numbers 



\ 2 / i ! \ 2 

n\ I n + k 



k \ k 



k=0 

was an essential tool in Apery's proof of the irrationality of 



CO 

«3> = E^ 



J 3 



> sumrecurs ion (binomial (n,k) ~2*binomial (n+k,k) ~2 ,k, A(n) ) ; 

(n + 2) 3 A(n + 2) - (3 + 2 n) (17 n 2 + 51 n + 39) A(n + 1) + (n + l) 3 A(n) = . 

Dougall's identity 



7-^6 



a, 1 + |, 6, c, d, 1 + 2a — 6 — c — d + n, —n 
|, 1 + a — 6, 1 + a — c, 1+a — d, a — n, 1+a+ra 

(1 + a) n (1 + a — b — c) n (1 + a — 6 — d) n (1 + a — c — d) n 
(1 + a - 6) n (l + a - c) n (1 + a - d) n (1 + a - 6 - c - d) n 

is proven by 

> sumrecursion(hyperterm( [a, l+a/2 ,b, c ,d, l+2*a-b-c-d+n, -n] , 

> [a/2 , 1+a-b , 1+a-c , 1+a-d , b+c+d-a-n , 1+a+n] , 1 , k) , k , S (n) ) ; 

— (a — d + n + l)(n + l + a — c) (n + 1 — 6 + a) (—b — c — d + a + n + 1) 
S(n + 1) + (1 + a + n) (n + 1 - c + a - d) (n + 1 - b + a - d) 
(n-c + o + 1-6) S(n) = . 
From this result, the right-hand side (|) of Dougall's identity can be read off 
directly. The complete computation is performed by 

> closedf orm(hyperterm( [a, l+a/2, b,c,d, l+2*a-b-c-d+n,-n] , 

> [a/2 , 1+a-b , 1+a-c , 1+a-d , b+c+d-a-n , 1+a+n] , 1 , k) , k , n) ; 

pochhammer(a + 1, n) pochhammer(— d — c + 1 + a, n) 

pochhammer(— d + 1 — b + a, n) pochhammer(a + 1 — b — c, n) /( 
pochhammer(l + a — d, n) pochhammer(l + a — c, n) 
pochhammer(l + a — b, n) pochhammer(a — d+1 — 6 — c, n)) . 
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Notice how important rational factorization is for such examples! 
The Wilson polynomials have the representation^ 

—n,a+b+c+d+n—l,a—x,a+x \ 
a + b, a + c, a + d J 

They include all classical systems like the Jacobi and Hahn polynomials. We get 

> sumrecursion(hyperterm( [-n,a+b+c+d+n-l ,a-x,a+x] , 

> [a+b, a+c, a+d] ,l,k) ,k,W(n)) ; 

(d + a + n + l)(n + l + a + c)(n + b + l + a)(a + 2n + c + b + d) 
(a + b + c + d + n) W(n + 2) - (2 n + 1 + a + b + c + d)(8 cdn 
+ 8bdn + bda 2 + 8bcn + 3b 2 n + bca 2 + 2b 2 d + Qanbc 
+ 7dn 2 + Qbncd + 3d 2 n + 2cd 2 + 2bd 2 + 7cn 2 + 3c 2 n 
+ 2c 2 d + 2bc 2 + 7bn 2 + 2b 2 c + 7 an 2 + 8adn + An 3 a 
+ 2d 2 n 2 + Adn 3 + 2c 2 n 2 + Acn 3 + 2b 2 n 2 + 46 n 3 + Qanbd 
+ b 2 + c 2 + 8acn + d 2 + 8abn + Aabc + 2ab + Q ancd 
+ bcd 2 + bc 2 d + ab 2 d + 2ab 2 n + acd 2 + 2ac 2 n + abd 2 
+ ab 2 c + abc 2 + b 2 cd + Aabcd + ac 2 d + 3 a 2 n + Qbn 2 a 
+ Qbn 2 d + Qbn 2 c + 2b 2 nd + 2b 2 nc + 4n 3 + 2 n 4 + 2c 2 nb 
+ 2 ad 2 + 2ac 2 + 2ab 2 + 6 cn 2 a + 6 cn 2 d + 2 c 2 n d + 2 a 2 cn 
+ 2a 2 bn + 2x 2 bd + 2x 2 cd + 2x 2 ad + x 2 a 2 + 2x 2 bc 
+ 2x 2 ac + x 2 d 2 + 2x 2 ab + x 2 b 2 + x 2 c 2 + Abed + 2d 2 nc 
+ 2d 2 nb + 2a 2 d + Qdn 2 a + 4abd + 2n 2 + 2a 2 b + 2a 2 c 
+ 2bc + 4acd + 2ad 2 n + a 2 + 2ac + 3bn + 2bd + 2cd 
+ 3dn + 2ad + 2a 2 dn + 3an + Aax 2 n + Abx 2 n + Acx 2 n 
+ 4dx 2 n + 4x 2 n 2 + 2ax 2 + 2bx 2 + 2cx 2 + 2dx 2 + 4x 2 n 
+ 3cn + 2a 2 n 2 + cd a 2 )W(n + 1) + (n + 1) (n + d + c) (n + b + d) 
(n + b + c)(d + 2n + a + b + c + 2) W(n) = , 

a recurrence equation for W n (x) which, however, is rather complicated since the 

middle coefficient admits no rational factorization. 

One knows from the theory that the recurrence equation has a special form 

which can be found by the command Sumrecursion: 

> Sumrecursion(hyperterm( [-n,n+a+b+c+d-l,a+x,a-x] , 

> [a+b, a+c, a+d] ,l,k) ,k,W(n,x)) ; 

(x — a) (a + x) W(n, x) = ((a + d + n) {n + a + c) {n + b + a) 
{a + b + c + d + n- I) W(n + 1, x)) /( 
(a + 2n + c + b + d)(2n + b + a + c + d-l))-( 
(a + d + n) (n + a + c) {n + b + a) (a + b + c + d + n — 1) 
(a + 2n + c + b + d) (2 n + b + a + c + d - 1) 
5 Sometimes a different standardization is used. But this is not essential. 
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WJx) 



4-^3 



+ 



(n + d + c-l)(n + b + d-l)(n + b-l + c)n 



(2n + b + a + c + d - 1) (2n - 2 + b + d + a + c) 
W(n, x) + 

(n + d + c-l)(n + b + d-l)(n + b-l + c)n W(n - 1, x) 

(2n + b + a + c + d - 1) (2n - 2 + b + d + a + c) ' 
A similar recurrence equation exists w.r.t. x: 

> Sumrecursion(hyperterm( [-n,n+a+b+c+d-l ,a+x,a-x] , 

> [a+b , a+c , a+d] , l,k) ,k,W(x,n)) ; 

(a + b + c + d + n — l)nW(x, n) = 

1 (x + d) (x + c) (x + 6) (a + x) W(x + 1, n) 

2 (2x + l)x 

1 (x + d) (x + c) (x + b) (a + x) 1 (x — d) (x — c) (x — b) (x — a) 



(2x + l)x 



+ 



'-l + 2x)x 



TTT . . 1 (x — d) (x — c) (x — b) (x — a) W(x — 1, n) 
W(x, n) + - — J -±— — — 

v ' ' 2 (-l + 2x)x 



Clausen's formula 

2^1 1 



a, b 
a+6+1/2 



x 



3 r 2 



2a, 2b, a + b 
a+b+l/2,2a + 2b 



x 



gives the cases when the square of a 2-^1 function is a 3_F 2 . The right-hand side 
is deduced from the left-hand side by 

> sumrecursion(hyperterm( [a,b] , [a+b+1/2] ,x,j)* 

> hyperterm( [a,b] , [a+b+1/2] ,x,k-j) , j ,C(k)) ; 

-(k + 1) (2 a + 1 + 2 b + 2 k) (2 a + 2 b + k) C(k + 1) 
+ 2 x (k + 2 b) (k + 2 a) (a + b + k) C(k) = , 
computing the coefficient of the Cauchy product. The resulting hypergeometric 
term can be obtained in one step byQ 

> Closedf orm(hyperterm( [a,b] , [a+b+1/2] , x, j)* 

> hyperterm( [a,b] , [a+b+1/2] ,x,k-j ) , j ,k) ; 

Hyperterm([2 6, 2 a, b + a], [a + b + —, 2 a + 2 b], x, k) 
The computation of a specific Feynman diagram yields the representation 



V(a,/9, 7 ) = (-l)° + ^ ■ T{a t P tZ 



■ d/2 ) r (d/2 - 7) r (a + 7 - d/2 ) r (0 + 7 - d/2 ) 



T(a)T((3)T(d/2)T(a + (3 + 2 7 - d)(m 2 ) a +P+"i- d 



a + (3 + 7 — d , a + r y — d/2 
a + {3 + 27 - d 



6 The Closedf orm procedure differs from the closedf orm procedure in that the hyperge- 
omtric term is not evaluated. 



22 



Since one is interested to compute this function for a, (3, 7 G N, and since the 
computation is easy for a, (3, 7 e {0, 1}, recurrence equations w.r.t. these variables 
can be used. Here is one w.r.t. (3: 

> sumrecursion( (-1) " ( alpha+bet a+gamma)* 

> GAMMA ( alpha+bet a+gamma-d/ 2 ) * GAMMA (d/2-gamma) * 

> GAMMA ( alpha+gamma-d/ 2 ) *GAMMA (bet a+gamma-d/ 2 ) / 

> (GAMMA (alpha) * GAMMA (beta) * GAMMA (d/2) * 

> GAMMA (alpha+beta+2*gamma-d) *(m~2) ~ 

> ( alpha+bet a+gamma-d) ) * 

> hyperterm( [alpha+beta+gamma-d, 

> alpha+gamma-d/2] , [alpha+beta+2*gamma-d] ,z,k) ,k,V(beta)) ; 

8 (3m 4 (13 + 1) (a + (3 + 1 + 7 - d) z V((3 + 2) + 2(3m 2 

(4-f + 2zp + 2z-zd + 2a + 2p-2d)(2a + 2(3 + 2 + 2-f-d) 
V(/3 + l) + 

(2 a + 2 (3 + 2-f - d) (2-f - d + 2 (3) (2 a + 2 (3 + 2 + 2-f - d) V((3) 
= . 

Similarly, one obtains recurrence equations w.r.t. a and 7. 

In some instances, Zeilberger's algorithm does not find the recurrence equation 
of lowest order. Assume, e.g., we want to deduce the right-hand side from the 
left-hand side of the identity 




Therefore, by Zeilberger's algorithm we compute a recurrence equation 

> RE:=sumrecursion( (-1) ~k*binomial (n,k) *binomial(3*k,n) ,k, S(n) ) ; 
RE :=2(3 + 2 n) S(n + 2) + 3 (7 + 5 n) S(n + 1) + 9 (n + 1) S(n) = 

and apply Petkovsek's algorithm to find its hypergeometric term solutions 

> rechyper(RE,S(n)) ; 

{-3} 

which gives the term ratio S n+1 /S n of the resulting hypergeometric term S n = 
("3) n . 

Here is another application of Petkovsek's algorithm: 

> RE : = (n+4) *s (n+2) +s (n+1 ) - (n+1) * s (n) =0 ; 

RE := (n + 4) s(n + 2) + s(n + 1) - (n + 1) s(n) = 

> rechyper (RE, s(n) ) ; 

f n + 1 (5 + 2n)(n+l) 

n + 3' (3 + 2n)(n + 3)^ ' 
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If continuous variables are involved, one can also compute holonomic differential 
equations for sums by a Zeilberger type algorithm which is implemented in the 
sumdif f eq procedure. 

We take some of the series representations of the Legendre polynomials to 
deduce the corresponding differential equation: 

> sumdif feq (binomial (n,k) *binomial (-n-1 ,k) * ( (l-x)/2) ~k,k,P(x) ) ; 

-(-1 + x) (x + 1) (i% P(x)) -2x{^- P(x)) + P(a;) n (n + 1) = 
ox z ox 

> sumdif f eq(l/2~n*binomial (n,k) ~2* (x-1) ~ (n-k) * (x+1) "k, k,P (x) ) ; 

-(-1+x) (x + 1) (-^-P(x))-2x(-^-P(x)) + P(x)n(n + l) = . 
ox ox 
To show a quadratic transformation like 



a, b 
2b 



Ax 



= (l + x) 2a - 2 F 1 



1 + x) 



a, a — b + 1/2 
6+1/2 



x 



we prove that both sides satisfy the same differential equation, and show that 
enough initial values agree. 

> sumdif feq(hyperterm( [a, b] , [2*b] ,4*x/(l+x)"2,k) ,k,Q(x)) ; 

-x (x - 1) (1 + xf (— Q(x)) + 2 (1 + x) {-x 2 + bx 2 - 2x a + b) (— Q(x)) 
ox z ox 

+ 4Q{x) {x-l)ab = 

> sumdiffeq((l+x)"(2*a)*hyperterm([a,a-b+l/2] , [b+1/2] ,x~2,k) ,k,Q(x)) ; 

-x (x - 1) (1 + x) 2 ( Q(x)) +2 (1 + x) (-x 2 + bx 2 - 2x a + b) Q(x)) 

+ AQ(x) (x-l)ab = 

> eval( [hypergeom( [a,b] , [2*b] ,4*x/(l+x) ~2)= 

> (l+x)~(2*a)*nypergeom([a,a-b+l/2] , [b+1/2] ,x"2) , 

> diff (hypergeom( [a,b] , [2*b] ,4*x/(l+x) "2)= 

> (l+x)~(2*a)*hypergeom([a,a-b+l/2] , [b+1/2] ,x~2) ,x)] ,x=0) ; 

[1 = 1, 2 a = 2 a] . 
On p. 258 in Ramanujan's second notebook one finds the identity 



2^1 



1 2 

3' 3 



X 



With Garvan we can ask the question: For which A, B, C, a, b, c, d is 



A,B 
C 



1 - 



1-x 
l + 2x 



[l + 2x) d 2 F 1 \ 



a, b 



,3 ? 



A computation with Maple gives 
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first :=hyperterm( [A,B] , [C] , l-( (l-x)/(l+2*x) ) ~3,k) ; 



pochhammer(A, k) pochhammer(.B, k) (1 



(1 - x 



1 3 



\k 



first := ( 1 + 2x r 

pochhammer(C, k) k\ 

> second:=(2*x+l) ~d*hyperterm( [a,b] , [c] ,x~3,k) ; 

(1 + 2x) d pochhammer(a, k) pochhammer(6, k) (x 3 ) k 

second := — — — 

pocnnammer(c, k) k\ 

> DEl:=sumdiffeq(first,k,S(x)) ; 

d 2 

DEI := x (x - 1) (1 + x + x 2 ) (1 + 2 x) 2 (— S(x)) + (1 + 2 x)(4x 4 + 9 B x 3 

ox 2 

+ 9Ax 3 -8Cx 3 + 3x 3 - 12 C x 2 + 9B x 2 + 9 Ax 2 + 3 x 2 -x 
- 6C x + 9 Ax + 9 B x - C)(^-S{x)) + 9 (x - l) 2 B AS{x) = 

> DE2:=sumdiffeq(second,k,S(x)) ; 



d 2 

DE2 := x (x - 1) (1 + x + x 2 ) (1 + 2xf (— S(x)) + (1 + 2x)(2a; 4 + 6bx 4 
-4dx 4 + Qax' k + 3bx 3 + 3ax 3 + x 3 -6cx + 4:X + 4:dx + 2 

Q 

-3c)(— S(x)) + (-12 x 4 ad + 36bax 4 + 4d 2 x 4 - Y2x 4 bd 
ox 

-2 x 3 d-6 x 3 ad + 36bax 3 -6 x 3 bd + 9bax 2 - 12 dx-Ad 2 x 
+ 12cxd + 6cd-4d)S(x) = 
DE:=collect (collect (lhs (DEI) -lhs(DE2) ,S(x) ) ,dif f (S(x) ,x)) ; 

DE := ((I + 2x)(Ax 4 + 9 B x 3 + 9 Ax 3 - 8C x 3 + 3x 3 - 12C x 2 + 9 B x 2 

+ 9 Ax 2 + 3 x 2 - x - QC x + 9 Ax + 9 B x - C) - (1 + 2 x)(2 x 4 

+ Qbx 4 -Adx 4 + Qax 4 + 3bx 3 + 3ax 3 + x 3 -Qcx + 4x + 4:dx 
d 

+ 2-3c))(— S(x)) + (9 (x- 1) 2 BA+ 12 x 4 ad - 36bax 4 
ox 

-4d 2 x 4 + 12x 4 bd + 2x 3 d + 6x 3 ad-36bax 3 + 6x 3 bd 
-9b ax 2 + I2dx + 4d 2 x- 12 cx d - 6 cd + 4 d)S(x) 
first coeff :=collect (f rontend(coef f , [DE, S(x)] ) ,x) ; 

firstcoeff := (12 a d - 36 ba - 4 d 2 + 12 bd) x 4 

+ (2d + 6ad-3Qba + Qbd)x 3 + (-9ba + 9B A)x 2 
+ (-18 BA+ \2d + Ad 2 - 12 cd) x + 9 B A - 6 cd + Ad 
secondcoef f : =collect (front end (coeff , [DE,dif f (S(x) ,x)] ) ,x) ; 

secondcoeff := (4 - 12 b + 8 d - 12 a) x 5 

+ (6 + 18 B + 18 A - 16 C + 4 d - 12 a - 12 b) x 4 
+ (-36- 3(1 + 8 + 275 + 27,4- 32 C) x 3 
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> 



+ (12 c- 7- 8d- 24 6 + 275 + 27 A) x 2 
+ (-9-8C + 9A + 9B-4d+12c)x-C-2 + 3c 
LIST:={coeff s(f irstcoeff ,x)} union {coeff s (secondcoef f ,x)}; 

LIST := {12ad-36ba- 4d 2 + 12bd, 9BA-6cd + 4d, 

-9ba + 9BA, -18 B A + 12 d + 4d 2 - 12 cd, 
6 + 18 B + 18 A - 16 6 + 4 d - 12 a - 12 b, 4 - 12 b + 8 d - 12 a, 
-36-3a + 8 + 275 + 27A-32 6, 
-9-8C + 9A + 9B-4d+12c, 
12c-7-8d-24 6 + 275 + 27A, -6-2 + 3c, 
2d + 6ad-36ba + 6bd} 
solve(LIST,{A,B,C,a,b,c,d}) ; 

1 



{d = d, B = ^d, A=^d+^,b = ^d, C =~ + ^d, a = ~d+~, c=^ + ~d}, 

' ' " . 1 1 , 1 j 

6 = — I — a, a = - a, 

9 ' 3 ' 



d Odd Z, 

1 1111 

\d = d, B = - d. A = - d+ -. b = -d + -. 

1 ' 3 3 3 3 3' 



c = 



6 

5 

C= 6 



2 2 

1 1 j 1 j 1 
- + -d, a = - a + 



■ d"}, {d = d, B = \ d + \, b = - d, C 

i 6 6 6 ^ ^ 

+ -d, A= -d},{d = d, B = -d + -,b = -d + -,C = - 



3' 

- d, 
2 ' 



1 j 5 1 . 1 , 

a =o d i c =a + a d i A= o d S 
6 6 6 3 

This leads to the unique solution 



'd/3,(l + d)/3 
(1 + 0/2 




since the hypergeometric functions are symmetric w.r.t. their upper parameters. 
Did you notice that this is exactly the computation from § 13/ To solve this 
question, a nonlinear system had to be solved. 

There is a theory of basic hypergeometric (g-hypergeometric) terms A\. for 
which Ak+i/Ak is rational w.r.t. q k . For almost all the results and algorithms 
corresponding q- versions exist. 

The corresponding series is called the basic hypergeometric series 



1 h b 2 



a, 
Ik 



q,x 



oo 

Y,A k x k 

k=0 



£ 

fc=0 



{a 1 ;q) k ■ (a 2 ; q) 



l)k---(a r ;q)kX k f k (*)W S 
■■■(K:a) k (a;a), V ; q J 



(bi;q)k(b2]q)k---(b s ;q)k(q; q) 
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where (a; q)k '■= Yi (1 — a g 7 ) is the q-Pochhammer-Symbol. is a q-hypergeometric 

j=0 

term and fulfils the recurrence equation {k e N) 



^4fc+i 



(l-aig fe )---(l-g r g fe ) 

(l-6 ig fc)...(l-6 8 gfc)(l-g*+l) 



•A 



with the initial value A := 1. 

All classical orthogonal polynomial families have at least one, most families 
possess several g-analogues. 

By the g-analogue of Zeilberger's algorithm, we get e.g. for the q-Laguerre 
polynomials 



4^) = ^^ 

(95 9)n 



q 



ct+l 



q, —xq 



n+a+1 



) 



the recurrence equation 

> read 'qsum.mpl' ; 

Copyright 1998, Harald Boeing & Wolfram Koepf 
Konrad — Zuse — Zentrum Berlin 

> qsumrecursion(qpochhammer (q~ (alpha+1) ,q,n)/qpochhammer(q,q,n)* 

> qphihyperterm( [q~ (-n)] , [q~ (alpha+1)] , q,-x*q~(n+alpha+l) ,k) , 

> q,k,L(n)); 

-q (-1 + q n ) L(n) + {-q 2 + g( 2n+a > x-q + g("+ a+1 ) + g( 1+ ™)) L(n - 1) 
+ q{q-q (a+n) )L(n-2) = . 
The g-analogue of Zeilberger's algorithm generates a third order recurrence equa- 
tion for the left-hand side of Jackson's g-analogue of Dixon's identity 



E 



n + b 




n + c 




~b + c~ 


n + k 


q 


c + k 


q 


b + k 



fc(3fc-l) 

q 2 



(?;?) 



n+b+c 



(?;?)„ (?;?) 6 (i'^)c 



> term : = (-1) ~k*qbinomial (n+b , n+k , q) *qbinomial (n+c , c+k , q) * 

> qbinomial(b+c,b+k,q)*q~(k*(3*k-l)/2) ; 

term := (— l) k qbinomial(n + 6, n + k, q) qbinomial(n + c, c + k, q) 
qbinomial(6 + c, b + k, q) q 0-/ 2k ( 3k - 1 )) 

> RE:=qsumrecursion(term,q,k,S(n)) ; 

RE : = -(-g (2n) + q) {q n + 1) (-1 + q n ) q 3 S(n) - (g 5 - g( 4 +"+ c + b ) + g 4 

q (3+2 n+b) _|_ ^(3+2 n+c) _ ^(3+2 n) _ q (3+c+b+n) 



q (3+n+b) _ q (3+n+c) 



q 3 + q (2n+b+2) + q (2n+c+2) 



_ q (2n+c+2+b) 



q (3n+c+2) _ q (3n+2) _ q (3n+b+2) 



_|_ q (b+c+2+4n) _ q (3n+l) + ? (l+4n+c+&) + ? (4n+c+6)^ 
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S(n - 1) + (g( 2+2 ") - g(3«+^) + q (2n + b+2) + g 4 



^(4+n+c+b) _ ^(3n+l+e+6) 



g (3+n+6) _|_ g (3+2n+6) + q 6 



_|_ g 5 _ g (3n+c+2+f>) _ + (2n+c+2) _ (3+n+c) 

+ g (3+2„+ C ) _ g (4+n+c))^ _ g (n+c+6)) _ 2 ) _ 

(g - g(" +c+6 )) (g 2 - q^) (q 2 - q^) (g 2 - q^+V) 
S(n-3) = . 

The g- analogue of Petkovsek's algorithm decides whether a g-holonomic recur- 
rence equation has g-hypergeometric term solutions. 

It finds the right-hand side of the g-analogue of Dixon's identity: 

> qrecsolve(RE,q,S(n) ) ; 

[[(g (n+1) - 1) S(n + 1) + (-gd+ot-Hn) + i) s(n) = 0]] . 
In many cases, much simpler is Paule's creative symmetrizing. With this method, 
we symmetrize the summand, and get the result in one step. 

> M : =qsimpcomb (subs (k=-k, term) /term, assume= [k, integer] ) ; 

M := q k 

> qsumrecursion((l+M)/2*term,q,k,S(n)) ; 

(1 - q n ) S(n) + (g(" +c+b ) - 1) S(n - 1) = . 
There exist similar algorithms for definite integration instead of summation. 
The Bateman integral representation 



4-C-l 



l-t) d -\F 1 



( a,b 






tx\ 


V c 





r( c + d) 



a, b 
c + d 



x 



is proven by 

> intrecursion(t~(c-l)*(l-t) "(d-l)* 

> hyperterm( [a,b] , [c] ,t*x,k) ,t ,B(k) ) ; 

-(k + l){k + d + c) B(k + 1) + B(k) x (b + k) (a + k) = 

> assume (d>0, c>0) ; 

> init:=int(t~(c-l)*(i-t)~(d-l) ,t=0. .1) ; 

init := B(c~ , d~) . 



We give our last example: On top of the g-Askey- Wilson scheme |13| we have the 
Askey-Wilson polynomials 

p n (x; a, 6, c, d\q) = (ab; q) n (ac; q) n (ad] q) n a~ n - 



4<P3 



q n ,abcdq n ,ae ,ae 



-i6 



ab, a c, a d 



q,Q 



The connection between those families can be written as 

n 

p n (x;a,(3,^,d\q) = ^ C m (n) p m (x;a,b, c, d\ q) 



m=0 
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with connection coefficients C m (n). 

Askey and Wilson showed the following representation for C m (n) 

c , n s = (®d,f3d,-fd,q;q) n (a /3 7 d q n ~ l ; q) m ^_ nm ^ 
(ad,(3d,jd,q,abcdq m - l ;q) m (q;q) n _ m 



q m ~ n , aP'jd q n + m ~\ adq m ,bdq m ,cdq r 
: '"' ] 1 abcdq 2m ,adq m ,(3dq m ^dq m 



In the general case this is not a g-hypergeometric term. 

However, for the special case (3 = b and 7 = c we get by the g-Zeilberger 
algorithm compare J7|, Sections 7.5 and 7.6) 

( y a/a;q) n _ m (abcdq n - 1 ;q) m ( y q,bc,bd,cd;q) n a n - m 

C m (n) 



(q, be, bd, cd; q) m {abcdq m ~ l ; q) m (q, abcdq 2m ; q) n _ m 

> c: = 'c': d:='d': 

> term : =qpochhammer (alpha*d , beta*d , gamma*d , q, q, n) * 

> qpochhammer (alpha*beta*gamma*d*q~ (n-1) ,q,m)/ 

> qpochhammer (alpha*d, beta*d, gamma*d, q, a*b*c*d*q~ (m-1) ,q,m)/ 

> qpochhammer (q,q,n-m) *q~ (m~2-n*m) *d~ (m-n) * 

> qphihyperterm( [q~ (m-n) , alpha*beta*gamma*d*q~ (n+m-1) , 

> a*d*q~m,b*d*q~m, c*d*q~m] , 

> [a*b*c*d*q~ (2*m) , alpha*d*q~m,beta*d*q~m,gamma*d*q~m] ,q,q,j): 

> qsumrecursion(subs({beta=b,gamma=c } ,term) ,q, j ,C(m) ) ; 

-(-1 + q m ) {-q + cdq m ) (caq {n+m) bd-q)(bcq m - q) {-q + bdq m ) 
(-q 3 + abcdq (2m) )(aq m -a q n ) C(m) + {-q 2 + caq m b d) 
{-q + abcdq {2m) )q 2 (-g m + q {1+n) ) (-g 2 + acq {n+m) bd) 
C(m-l) = 

and similar results for a = a, 7 = c and for a = a, (3 = b: 

From these results one can derive the connection coefficients between many 
families of the g-Askey- Wilson tableau by limit computations. 
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